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We review the basic theory of matrix product states (MPS) as a nu- 

Smerical variational ansatz for time evolution, and present two methods 
to simulate finite temperature systems with MPS: the ancilla method 
and the minimally entangled typical thermal state method. A sample 
C"j calculation with the Bose-Hubbard model is provided. 

o 

1.1. Introduction 

> 

pTj The dimension of the Hilbert space for a general many-body system in- 

creases exponentially with the system size, severely restricting the sizes 
which are amenable to straightforward numerical study. Several techniques 
have been developed to deal with this fact, such as the stochastic sam- 
pling of the Hilbert space in quantum Monte Carlo and the judicious use 



00 

o 

of symmetries and sparse matrix structures in exact diagonalizations. The 
most successful approximate method for Id systems is the density matrix 
renormalization group (DMRG) method first pioneered by White (T). Soon 
after, the theory of matrix product states [2l [3] (MPS) was used to shed 
light on the amazing success of DMRG [4j |5j. Ideas from quantum infor- 
mation theory, most notably the idea of bipartite entanglement, have led to 
the development of MPS algorithms which generalize DMRG to time evolu- 
tion [6j [7] , periodic boundary conditions [8] , and finite temperature [9 10 



A thorough discussion of the time-evolving block decimation algorithm, an 
MPS algorithm for zero temperature time evolution, is given in chapter 
??. In this chapter we review algorithms based on MPS for finite tempera- 
ture simulations and discuss their relevance to studying finite temperature 
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superfluid systems. 

1.2. Methodology 

1.2.1. Matrix Product States 

A matrix product stat^ (MPS) is defined as 



l^mps) 



»i,i2,...i£=l 



TrfA^-.-A^) •••,»£> 



(1.1) 



where the A' k ' lk are matrices^ the dimension of which is a fixed number 
X known as the bond dimension, d is the dimension of the Hilbert space 
spanned by the {|*fc)}, and L is the number of lattice sites. Let us refer 
to the set of all MPSs with bond dimension \ as An MPS in A4 X 

contains Ld\ 2 parameters, and so it is clear that any state on a finite 
lattice can be written as an MPS provided we take the bond dimension to 
be Xmax = rfL L / 2 J . However, the great utility of MPSs is that an MPS with 
bond dimension \ ^ Xmax often provides an excellent approximation to 
the true state [12] and allows for exponentially more efficient manipulation 
and calculation of observables than an exact representation. 
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Fig. 1.1. (a) Tensor network representation of full 4 site wavefunction. (b) Tensor 
network representation of an MPS on 4 sites, (c) Tensor network representation of an 
MPO on 4 sites. 



To visualize MPSs and operations with them, it is useful to introduce 
the notion of a tensor network diagram as in Fig. |1.1| In such a diagram a 
box represents a tensor, free lines are uncontracted indices and closed lines 

a An MPS is a vector in Hilbert space. The qualifier matrix product refers to the fact 
that the expansion coefficients in the Fock basis are expressed as products of matrices. 
b These matrices can be taken to have the same symmetry as the state they represent, 
e.g., if the state has real coefficients in some basis then the MPS matrices can be taken 
to be real. See [11| and references therein for more details. 
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are contracted indices. Fig. |l.l| a) shows the state of a many-body system 
expressed in the basis of the full Hilbert space as an L-index tensor, and 
Fig. shows the same state written as an MPS. The advantage of the 

MPS representation becomes clear when we compute scalar products such 
as (1>\d\4>). 

Before we discuss how scalar products are efficiently computed, it is 
advantageous to introduce a matrix product operator (MPO) as 



6= J2 E ^(MW^-.-MW^Iii.-.-.iiXii, 

ii,...,iz,=l i[,...,i' L = l 



■■■ ,*L\ 



where each of the M[ k l lkl k is a matrix the dimensions of which are bounded 
by a fixed number xo known as the bond dimension. The tensor network 
representation of an MPO is similar to that of an MPS, but there are two 
uncontracted indices per tensor corresponding to the bra and ket indices; 
see Fig. |l.l[ c) . Equivalently, one can think of each element of the matrix 
as being operator valued, where the operator acts on the space spanned 
by{|«fc)}- 

Let us now see how to evaluate the scalar product of an operator O 
represented as an MPO between two states and \<f>) represented as MPSs. 
Let us denote the MPO matrices of O as M and the MPS matrices of \ip) 
and \(f>) as A and B, respectively. Then, we have 



101. 



]T E ^(aW^.-aM^Ti^mW^ 

ii,...,it=l i[ , . . . ,i' L = l 



xTr(B [1]il ...B [L]iL 

d 



(1.2) 



Tr 



A [1]il * ® M [1Iili 'i ® B [1]i 'i 

d 

A [LIiL *(g)M [L]iLi L(g)B [Lli L 



Tr 



A, B) 



• eW(a,b 



(1.3) 
(1.4) 



where the last line defines the generalized transfer matrix E^j (A, B) = 
Ya i'=i AM ik *(g)MM iki k(8BM i k, which is a \ 2 Xo x x 2 Xo matrix. Naively 
we would expect that the multiplication of two transfer matrices would 
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require O (x 6 Xo) operations, but the special structure of the transfer ma- 
trices allows us to perform such a multiplication in O (x 5 Xo c ' 2 )|2] as 



e« (a,b)e& +1] (a,b; 



d XO 



EE EE ([ g m(a,b 



i' = l/9=l \i=l 7 "=l 



G*m (A,B) 



a7/3],[a'7"/3'] 



M [k+l]ii' 

7"7' 



B 



(1.5) 

[k+i]i' 



[a7/3],[a'7"/3'] 



E [eS(a-b 



[a7/3],[""P7"/3"] 



Here the square brackets around indices denote a composite index in the 
Kronecker representation and parentheses give the order in which the con- 
traction should be performed to ensure the best scaling. In particular, it 
is essential not to sum over the a" and j3" indices simultaneously]^] The 
tensor network representation of the scalar product procedure is given in 
Fig.Ol 
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Fig. 1.2. Tensor network representation of the scalar product procedure Eqs. | |l,3[ l- |l~4"| l. 
The transfer matrices (A, B) have been abbreviated as for succinctness. 



Many operators of interest, such as translation invariant Id Hamiltoni- 
ans, can be easily represented as MPOs with small bond dimension \o ~4- 



10 15 16 , and the MPO representations of more complex operators can 



c The fact that the boundary matrices of MPSs with open boundary conditions have 
bond dimension 1 allows us to perform this contraction in 0(x 3 Xo^ 2 )> an< ^ recent 
developments for periodic boundary conditions have reduced the scaling to O (x 3 Xo^ 2 ) 
for large systems with only a few relevant correlation lengths [131 114] . 
d Here and throughout we use greek indices to denote bond indices and roman indices to 
denote physical indices. 
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be constructed using simple MPO arithmetic 15 17 . That the MPO form 



of an operator is optimal for MPS algorithms can be straightforwardly de- 
duced using the tensor network formalism, as the scalar product of an MPO 
between two MPSs is the most general tensor network that can be efficiently 
contracted; see Fig. \1.2\ 

We now turn to the simulation of time evolution using MPSs. The 
main difficulty of using MPSs is that M. x is not a vector spaceQ Thus, 
when operators such as the propagator are applied to an MPS we must hnd 
the optimaj^] projection into M. x to keep the algorithm efficient. We denote 



this projection as V x . The optimal MPS 

U\<f>) is 



G M. x representing the MPS 



u\4>) 



mm 

min 

|V>eA<» 



\1>)-U\<f>) 

\^) + (<t>\tfu\ ( l>)-m(^\u\<t>) 



(1.6) 
(1.7) 



where 1Z(») denotes the real part. Each of the scalar products in Eq. (1.7) 



may be written as a quadratic form in each of the matrices Al k l Ik , as is 
demonstrated in the tensor network diagram Fig. |1.3[ 




Fig. 1.3. Tensor network representation of the quadratic form representing (ip\U\4>) in 
Eq. ( fT~7| l. 



e This can be seen from the fact that the addition of two MPSs is given by the direct 
sum of their matrices: |V>c) = \4>a) + IV'b) => C' k ' = AW ffi BW. If the matrices A" 
and BW have orthogonal bases then dim (C' k l) = dim (AM) + dim (BW). 
f By optimal we mean that the overlap is maximal in the 2-norm. Although MPSs do 
not form a vector space, they are embedded in a larger Hilbert space and so this norm 
is well-defined. 
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Again denoting the matrices in the MPS representation of by A and 
those of |0) by B, the quadratic form of the k th site may be written as 

= A^QMAM + B^qM -BM - 2TZ (a^Q^bM) , (1.8) 

where A^l represents the d% 2 elements of the |AM lk }, arranged as a vector, 
and the matrices Qq are defined as 



Q 



XO 



[ai k a'] [/3iJ.i3'j 



7,7' = 1 



S(C,D) 



(1.9) 



[c<7/3],[a'7'/3'] 

where C and D are either A or B depending on the quadratic form. The 
M Ik ' k in this final expression are the matrices in the MPO representation 
of O. The stationary points of the quadratic form Eq. (1.8) are given by 
the solution of the linear systenjs] 1 8 

QjAM = Q^bM . (1.10) 

The algorithmic procedure for time evolution is to sweep back and forth 



through the lattice, solving Eq. (1.10) at each site until convergence is 
reached. In practice, it is essential for efficiency not to explicitly form 
the matrices Q # , but rather to use iterative methods which require only 
multiplication by the Q, to solve Eq. (1.10). Details on the form of the 



propagator U can be found in 17 18 



1.2.2. The Ancilla Method 

At finite temperature, the state of a quantum system is g iven by the thermal 
density matrix p = e~@ H /Z. The ancilla method |9l|l9j relies on the notion 



of purification 20 to represent the thermal density matrix as a pure state in 



an enlarged Hilbert space. Each physical site is augmented with an ancilla 
which has the same Hilbert space dimension as the physical site. The MPS 
representation of such a system is 



d 

E 

i 1 ,...,i L = l Oi, 



d 

E 



Tr A 



. A[ L l iLaL 



\ixai ■ ■ ■ ilQ>l) 



(1.11) 



g It is important to note that while Qj is the quadratic form representing the scalar 
product (ip\ip) it can not in general be made equal to the identity. The numerical 
conditioning of this matrix and of the linear system Eq. (flTTol can be improved by 
suitable choice of "gauge conditions" on the matrices A; see [81. 
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One can think of the combined system as a two-legged ladder, with the 
physical sites on the lower leg and the ancillae on the upper leg. The pur- 
pose of the ancillae is to act as a perfect heat bath which, when traced out, 
provides the proper thermal density matrix for the physical system. The 
choice of ancilla for infinite temperature (/3 = 0) is simply the normalized 
purification of the identity 

IV> (0)) = ^ nti E£,«»=l Si k a„ %a k ) , (1.12) 

which represents a product of maximally entangled site-ancilla pairs. This 
state has an MPS representation with bond dimension 1 generated by tak- 
ing all matrices to be A^L lkak = <5 Qi i^ j i<5i kak /v / d. The extension to finite 
inverse temperature (3 is provided by evolving only the physical sitesj^] in 
imaginary time up to /3/2, 

hM/3)>=e-^/ 2 KM0)). (1-13) 
This time evolution can be efficiently performed using the methods of 



Sec. 1.2.1 Observables are calculated using transfer matrices as above with 
the additional requirement that the ancilla degrees of freedom are traced 
over. 

The ancilla method is conceptually very simple, and becomes numeri- 
cally exact for large enough bond dimension. However, because the MPS 



Eq. ( 1.11 ) has to encode the information of both the system and the bath, 
it requires a bond dimension ~ 7Q. S . & t low temperatures, where Xg.s. is the 
bond dimension required to accurately represent the ground state. Typical 
values of Xg.s. range from 50-5000, making the ancilla method impractical 
for many systems at very low temperatures. 

We conclude this section by remarking that the ancilla method repre- 
sents a highly idealized heat bath chosen to reproduce the exact thermal 
density matrix. Many of the current examples of strongly correlated many- 
body systems, e. g. cold atoms, are very mesoscopic and are in contact with 
thermal reservoirs which are also mesoscopic. A modification of the ancilla 
method where the perfect entanglement at infinite temperature is replaced 
with ancilla-ancilla and ancilla-system couplings in the Hamiltonian can be 
devised. Alternatively, one can directly simulate master equations by con- 
sidering matrix product density operators with optimal projections based 



or matrix product decompositions of 



on the Hilbert-Schmidt distance 
superkets with local projections [10 

h That is, the Hamiltonian only couples physical sites to physical sites, and not ancillae 
to ancillae or physical sites to ancillae. 
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1.2.3. Minimally Entangled Typical Thermal States 

A new method for finite temperature MPS simulations has recently been 



proposed by White 21 . The idea stems from the question "What is a 
typical wave function of a quantum system at finite temperature?" That 
is, if we are to measure a quantum system at finite temperature, what 
"typical" pure states would we find, and with what probabilities? It is 
clear from the basic tenets of statistical mechanics that any set of typical 
states {\4>{i))} must satisfy 

J>(i)l0WW(i)l = «r^, (1.14) 

i 

where P (i) is the probability of measuring the system to be in state \<p(i)), 
and so the expectation of an operator A at finite temperature may be 
written as 

(i> = £^<<M*)|i|0«>, (i-iB) 



with Z the partition function. From Eq. (1.15), we see that we can calcu- 
late observables using an unweighted average of {<j> (i) \A\</> (i)) if we choose 
the \4> (i)} at random according to their probabilities of being measured, 
P (i) /Z. It is easy to generate states satisfying the typicality condition 
Eq. ( |1.14 1 simply by defining any orthonormal basis {\i)} and defining 



\<j> (<)) = [P (i)]- 1 ' 2 exp (-/3#/2) \i) ,P(i) = (i\ exp \i) . (1.16) 

We now use the freedom in the choice of the orthonormal basis {\i}} to 
generate typical states with the least amount of spatial entanglement, as 
these are the states which can be most efficiently represented as MPSs 
[6| [22]. This amounts to taking the to be classical product states 

(CPSs), |z) = rifc=i l*fc)i where labels the state of site k. The set of 
\4> (i)) obtained from this choice of {\i)} are called minimally entangled 
typical thermal states (METTS). 

The most efficient algorithmic procedure for generating thermal aver- 
ages using METTS is as follows: 

(1) Choose a CPS \i) at random. 



(2) Evolve in imaginary time using the methods of Sec 1.2.1 to generate 
the METTS \(f>(i)) = [P (i)P 1/2 exp (-0H/2) \i). 

(3) Compute observables of interest using this METTS and add to the 
running averages. 
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(4) Randomly select a new CPS \i') according to the probability \(i'\(f> (i))\ 2 - 

(5) Repeat from step 2 until converged. 

We see that the main loop of this algorithm closely resembles a Monte Carlo 
iteration with measurement taking the place of the usual configuration up- 
dates. However, it does not rely on a rejection method to perform sampling, 
and so each METTS that is generated can be used to generate statistics. 
In practice very few (~ 100) METTS suffice to obtain the total energy to 
a relative accuracy of 10 -5 . For algorithmic details on how to perform the 
CPS selection to minimize correlations between successive METTS we refer 
the reader to [17]. 

This METTS algorithm has many advantages over the ancilla method 
of the previous section. As we do not have to encode the bath degrees of 
freedom in our MPS, the bond dimension required to accurately represent 
each METTS ranges from 1 at infinite temperature to Xs-s. at very low 
temperatures. This makes the METTS method more efficient than the 
ancilla method by a factor of 10 3 — 10 10 for typical systems at very low 
temperatures. Additionally, if the Hamiltonian of interest has a global 
symmetry then we can use the fact that the MPS matrices must transform 
irreducibly to speed up the calculation [15] or find the thermal ensemble 
corresponding to a fixed quantum number (canonical ensemble). This latter 
point is relevant to cold atom systems where the total number of atoms is 
held fixedQ 



1.3. Validity issues 

It has been shown that MPSs can faithfully represent ground states of 
Id gapped Hamiltonians with at most nearest neighbor interactions with 



a bond dimension which grows only polynomially in the system size 12 



In higher dimensions this polynomial scaling gives way to an exponential 



scaling 24 , but calculations on 2D systems of width 8-12 are still feasible 



25 . Generalizations of MPSs to higher dimensions exist, but are so far 



limited by poor polynomial scaling of tensor contractions 26 -28 . Perhaps 



the most important quality of MPS methods as compared to other efficient 



'The ancilla method can also be used to simulate systems in the canonical ensemble, but 
the process is complicated by the fact that we need the purification of the constrained 
infinite temperature density matrix. This purification can be generated using a ground 
state DMRG-type calculation with a suitably chosen Hamiltonian [23] . The Hamiltonian 
will contain artificial ancilla-ancilla and ancilla-physical site couplings which are typically 
highly nonlocal. 
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many body methods, such as quantum Monte Carlo, is that MPS methods 
work equally well for fermionic or frustrated systems. All of the methods 
presented here will work equally well for any Id or quasi-ld physical system. 

1.4. Application: Specific Heat of the Hard-core Extended 
Bose-Hubbard Model 

As an example of how the above methods may be applied to study the 
behavior of a finite temperature supcrfluid system, we study the properties 
of the hard-core extended Bose-Hubbard model 

H = -JE {i , j) (S&+H.C.) +VZ {iJ) Wj (1-17) 

at half filling. This model is known to have a supcrfluid phase in the 
XY universality class for V < 2J. In the below figure we show a typical 
thermodynamic quantity, the specific heat Cy = f3 2 ((H 2 ) — (H) 2 )/L, as 
a function of temperature and the nearest-neighbor repulsion. Note that 
computation of (H 2 ) is easily performed when the MPO representation of 
H is known. 




Fig. 1.4. Specific heat of the hard-core extended Bose-Hubbard model on 34 sites for 
repulsive nearest neighbor interaction V = 0,0.5,1.0,1.5,2.0,3.0,4.0. The qualitative 
behavior of the low temperature specific heat changes as V becomes larger than 2J 
because the system transitions from a gapless superfluid phase into a gapped insulating 
phase. 
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